G4: Edinaldo de Alencar / Igor Freire / Ramon Araújo / Ricardo Ribeiro
29 de outubro de 2014
Pacote utilizado para leitura em formato de séries temporais:
library(zoo)
Importar as bases de dados utilizando read.zoo:
# A base de dados e os scripts R estão no mesmo diretório (o diretório atual)
setwd(paste("~/Documents/Mestrado/UFPA/Mineração de Dados/data-mining-ppgee/trabalho-2-forecast/", sep=""))
#setwd(paste("./", sep=""));
# Inicialmente, ler a base de dados diários como um data frame (através de read.csv)
dataframe_diario <- read.csv(file = "dataset_diario.csv", sep = ";", dec = ",", header = TRUE)
# Em seguida, converter para uma série temporal (lista indexada pela data)
dataset_diario <- zoo(as.matrix(dataframe_diario[, -1:-2]), as.Date(dataframe_diario[,1], format = "%d/%m/%y"))
# Obs: em "format" usa-se y minusculo, pois a data está no formato dd/mm/yy
# A função read.zoo() abaixo não retorna lista com vetores categóricos e numéricos ao mesmo tempo.
# Isto é, se houver dados categóricos e numéricos na base, os numéricos serão convertidos.
# Por isso, a base com dados diários não foi lida diretamente com read.zoo. A base com dados mensais
# pode ser lida diretamente com read.zoo()
# Ler o .csv como uma série temporal (indexado pela data)
dataset_mensal <- read.zoo(file = "dataset_mensal.csv", sep = ";", dec = ",", header = TRUE,
index = 1, tz = "", FUN = as.yearmon, format = "%m/%Y", drop = FALSE)
# index -> coluna do arquivo .csv que contém a data
# Obs: em "format" usa-se Y maiúsculo, pois a data está no formato mm/yyyy
dataset_diario é a base com valores de fluxos diários de 2002 a 2009.dataset_mensal é a base com os valores de fluxo mensais de 1992 a 2009.sep =)dec =)Observar que a informação sobre os dias da semana é redundante, pois pode ser obtida através de:
dia <- weekdays(index(dataset_diario))
que tem como resultado, por exemplo, para as 10 primeiras amostras da base:
head(dia, 10)
## [1] "Terça Feira" "Quarta Feira" "Quinta Feira" "Sexta Feira"
## [5] "Sábado" "Domingo" "Segunda Feira" "Terça Feira"
## [9] "Quarta Feira" "Quinta Feira"
o qual pode ser confirmado comparando com as 10 primeiras entradas na coluna “dia” da base:
head(dataframe_diario[,2], 10)
## [1] ter qua qui sex sab dom seg ter qua qui
## Levels: dom qua qui sab seg sex ter
Por isso, na leitura da base em dataset_diario, esta coluna foi ignorada.
summary(dataset_diario)
## Index dataset_diario
## Min. :2002-01-01 Min. :11130
## 1st Qu.:2004-01-01 1st Qu.:19468
## Median :2005-12-31 Median :22811
## Mean :2005-12-31 Mean :23158
## 3rd Qu.:2007-12-31 3rd Qu.:26712
## Max. :2009-12-31 Max. :34292
summary(dataset_mensal)
## Index fluxo
## Min. :1992 Min. :261544
## 1st Qu.:1996 1st Qu.:408961
## Median :2001 Median :529912
## Mean :2001 Mean :548440
## 3rd Qu.:2005 3rd Qu.:662250
## Max. :2010 Max. :982708
Conversão dos dados para o formato de série temporal no R.
# Frequency --> número de observações por unidade de tempo
# define a unidade de tempo (e.g. 12: unidade de tempo = ano)
tsMensal <- ts(dataset_mensal, frequency=12, start=1992)
plot(tsMensal, xlab="Anos", ylab="Fluxo Mensal")
Série temporal diária
tsDiario <- ts(dataset_diario, start=c(2002,1,1), frequency=365)
plot(tsDiario, xlab="Anos", ylab="Fluxo")
Decompor a série temporal em:
plot(decompose(tsMensal), xlab="Anos")
plot(decompose(tsDiario), xlab="Anos")
Realizar prospecções de curto, médio e longo prazo.
Para curto prazo, será utilizada a série temporal com dados diários (tsDiario).
Para médio e longo prazo, será utilizada a série temporal com dados mensais (tsMensal).
Separação das séries temporais em conjuntos de treinamento e de testes.
Separação de tsDiario em conjuntos de treinamento e teste.
tsDiarioTrain <- window(tsDiario, end=c(2005,366)) # De Jan 2002 a Dez 2005
plot(tsDiarioTrain, xlab="Anos", ylab="Fluxo Mensal")
Série temporal diária multi-sazonal
# Considerando sazonalidade semanal e anual
tsDiarioMSzl <- msts(tsDiario, seasonal.periods=c(7,365.25))
Série temporal multi-sazonal do conjunto de treinamento
# Considerando as sazonalidades semanal e anual
tsDiarioTrainMszl <- msts(tsDiarioTrain, seasonal.periods=c(7,365.25))
tsDiarioTest30Days <- window(tsDiario, start=c(2006,1), end=c(2006,30)) # De 1/1/2006 a 30/1/2006
# 45 dias a partir De 1/1/2006
tsDiarioTest45Days <- window(tsDiario, start=c(2006,1), end=c(2006,45))
tsMensalTrain <- window(tsMensal, end=c(2005, 12)) # De Jan 1992 a Dez 2005
Separação de tsMensal em conjuntos de treinamento e teste:
tsMensalTest4mth <- window(tsMensal, start=2006, end=c(2006,4)) # De Jan 2006 a Mar 2006
tsMensalTest6mth <- window(tsMensal, start=2006, end=c(2006,6)) # De Jan 2006 a Jun 2006
tsMensalTest1yr <- window(tsMensal, start=2006, end=2007) # De Jan 2006 a Jan 2007
tsMensalTest2yr <- window(tsMensal, start=2006, end=(2008)) # De Jan 2006 a Jan 2008
Pacote utilizado (ver (R. J. Hyndman 2014) e (Leek 2014)):
library(forecast)
Apresenta-se a seguir:
Nota:
Saída: objeto forecast, contendo:
Série Original
Predições
Método utilizado
Intervalos de Predição
Resíduos
Métricas comparativas usuais:
Diferentemente do ME, MAE e RMSE, o medidor MAPE é independente da escala.
\[M = \frac{1}{n} \sum \limits_{t=1}^{n} \left| \frac{A_t - F_t}{A_t} \right|,\]
onde \(A_t\) é o valor verdadeiro, \(F_t\) é o valor predito e \(n\) é o número de amostras temporais preditas.
Como desvantagem, o critério MAPE tem comportamento indefinido quando o valor verdadeiro \(A_t=0\) (divisão por zero).
Funções utilizadas na biblioteca forecast:
meanf()naive()snaive()rwf()Prospecções através dos métodos simplísticos para um prazo 30 dias.
Fluxo diário de Janeiro de 2002 a Dezembro de 2005
diasAPrever <- 30
# Média
f_mean <- meanf(tsDiarioTrain, h=diasAPrever)
# Naïve
f_naive <- naive(tsDiarioTrain, h=diasAPrever)
# Naïve Sazonal
f_seasonal_naive <- snaive(tsDiarioTrainMszl, h=diasAPrever)
# Drift
f_drift <- rwf(tsDiarioTrain, drift = TRUE, h=diasAPrever)
Acurácias MAPE obtidas.
## Mean Näive N-Sazonal Drift
## Training set 10.47 6.38 9.030 6.382
## Test set 12.66 10.15 7.664 9.951
## Mean Näive N-Sazonal Drift
## Training set 22.82 4.376 7.673 4.408
## Test set 29.64 5.202 6.735 6.139
## Mean Näive N-Sazonal Drift
## Training set 22.82 4.376 7.673 4.408
## Test set 34.45 5.524 7.163 4.601
ARIMA (Autoregressive Integrated Moving Average)
Combinação dos modelos:
Em um modelo de auto-regressão, prevemos a variável de interesse usando uma combinação linear dos valores passados da variável.
Assim, um modelo auto-regressivo de ordem \(p\), denotado por \(AR(p)\), pode ser escrito como:
\(y_t = c + \varphi_1 y_{t-1} + \varphi_2 y_{t-2} + \dots + \varphi_p y_{t-p} + e_t,\)
onde \(c\) é uma constante, \(e_t\) o ruído branco, e, os valores defasados de \(y_t\) são os preditores.
Para um modelo \(AR(1)\), tem-se \(\varphi_1\in(-1,1)\)
Para um modelo \(AR(2)\), tem-se \(\varphi_1,\varphi_2\in(-1,1)\) com \(\varphi_1 + \varphi_2 < 1\) e \(\varphi_2 - \varphi_1 < 1\)
Ao invés de usar valores passados da variável de previsão em uma regressão, o modelo de médias-móveis usa erros de previsão do passado em um modelo de regressão-like.
\(y_t = c + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} + \dots + \theta_q e_{t-q}.\)
Esse modelo tem ordem \(q\) e é denotado por \(MA(q)\).
Para um modelo \(MA(1)\), tem-se \(\theta_1\in(-1,1)\)
Para um modelo \(MA(2)\), tem-se \(\theta_1,\theta_2\in(-1,1)\) com \(\theta_1 + \theta_2 > -1\) e \(\theta_1 - \theta_2 < 1\)
É possível escrever o modelo \(AR(p)\) como um modelo \(MA(\infty)\). Em particular para um modelo \(AR(1)\), vem:
\(y_t = \varphi_1 y_{t-1} + e_t\)
\(y_t = \varphi_1 ( \varphi_1 y_{t-2} + e_{t-1} ) + e_t\)
\(y_t = \varphi_1^2 y_{t-2} + \varphi_1 e_{t-1} + e_t\)
\(y_t = \varphi_1^3 y_{t-3} + \varphi_1^2 e_{t-2} + \varphi_1 e_{t-1} + e_t\)
e assim por diante. Onde \(\varphi_1\in(-1,1)\)
mesesAPrever <- 4
Modelo:
arima_model_mensal <- auto.arima(tsMensalTrain)
Predição:
fcast_arima_4mth <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_4mth, xlab="Anos", ylab="Fluxo Mensal", include=36)
lines(tsMensalTest4mth, col="red")
Acurácia:
accuracy(fcast_arima_4mth, tsMensalTest4mth)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set 145 10236 7610 -0.0156 1.595 0.2052 -0.005073 NA
## Test set 2173 10234 8719 0.2861 1.276 0.2351 -0.511317 0.1839
mesesAPrever <- 6
Predição:
fcast_arima_6mth <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_6mth, xlab="Anos", ylab="Fluxo Mensal", include=36)
lines(tsMensalTest6mth, col="red")
Acurácia:
accuracy(fcast_arima_6mth, tsMensalTest6mth)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set 145 10236 7610 -0.0156 1.595 0.2052 -0.005073 NA
## Test set 3437 9902 8148 0.4719 1.181 0.2197 -0.352187 0.2011
mesesAPrever <- 12
Predição:
fcast_arima_1yr <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_1yr, include=48)
lines(tsMensalTest1yr, col="red")
Acurácia:
accuracy(fcast_arima_1yr, tsMensalTest1yr)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set 145 10236 7610 -0.0156 1.595 0.2052 -0.005073 NA
## Test set 22409 31001 24765 2.9237 3.278 0.6677 0.697780 0.7174
diasAPrever <- 30
Modelo, a partir da série multi-sazonal do conjunto de treinamento:
tbats_model <- tbats(tsDiarioTrainMszl)
Predição:
tbats_fc30Days <- forecast(tbats_model, h=diasAPrever)
plot(tbats_fc30Days, include = 120)
lines(tsDiarioTest30Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)
Acurácia:
accuracy(tbats_fc30Days, tsDiarioTest30Days)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set 4.494 687.8 393.3 -0.109 2.095 0.2101 0.0006724 NA
## Test set -210.910 754.8 600.6 -1.086 2.710 0.3209 0.7607279 0.4361
diasAPrever <- 45
Predição:
tbats_fc45Days <- forecast(tbats_model, h=diasAPrever)
plot(tbats_fc45Days, include = 120)
lines(tsDiarioTest45Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)
Acurácia:
accuracy(tbats_fc45Days, tsDiarioTest45Days)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set 4.494 687.8 393.3 -0.1090 2.095 0.2101 0.0006724 NA
## Test set -13.861 722.4 595.0 -0.2168 2.638 0.3179 0.7639143 0.3844
Ver (R. J. Hyndman and Athanasopoulos 2013). O autor possui um livro exclusivamente sobre o assunto em (R. Hyndman et al. 2008).
Para cada um dos 15 métodos, há dois modelos possíveis: com erro aditivo e com erro multiplicativo.
ets(ts) sem argumentos além da série temporal ts determina automaticamente o método apropriado.ets apresenta o modelo escolhido e o AIC resultante.ets não trabalha com séries cuja sazonalidade é superior a 24 unidades temporais. Portanto, será utilizada somente pra prospecções realizadas com dados mensais (sazonalidade = 12).ets parametrizado automáticamente deve dar resultados muito precisos para prospecções de poucos pontos (por exemplo, 4 pontos).mesesAPrever <- 6
Modelo:
etsMensal <- ets(tsMensalTrain)
Prospecção:
fcastMensal6mth <- forecast(etsMensal, h=mesesAPrever)
Comparar predição com valores reais do conjunto de teste:
plot(fcastMensal6mth, xlab="Anos", ylab="Fluxo Mensal", include = 36);
lines(tsMensalTest6mth, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)
Acurácia:
accuracy(fcastMensal6mth, tsMensalTest6mth)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set 910 8952 6764 0.1729 1.438 0.1824 -0.005476 NA
## Test set 4511 10353 10213 0.6550 1.495 0.2753 0.152408 0.1901
mesesAPrever <- 12
Prospecção:
fcastMensal1yr <- forecast(etsMensal, h=mesesAPrever)
Comparar predição com valores reais do conjunto de teste:
plot(fcastMensal1yr, xlab="Anos", ylab="Fluxo Mensal", include = 48);
lines(tsMensalTest1yr, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)
Acurácia:
accuracy(fcastMensal1yr, tsMensalTest1yr)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set 910 8952 6764 0.1729 1.438 0.1824 -0.005476 NA
## Test set 21648 29250 24499 2.8481 3.268 0.6605 0.789771 0.6718
Pode-se observar a existência de um fator multiplicativo que não está sendo capturado pelo modelo. Uma transformação Box-Cox pode ser apropriada.
Apropriado quando dados apresentam diferentes variações em diferentes níveis da série.
\(W_t = \begin{cases} \log\left(y_t\right) & \lambda = 0\\ \frac{(y_t^\lambda - 1)}{\lambda} & \lambda \neq 0\end{cases}\)
plot(tsMensalTrain, xlab="Anos", ylab="Fluxo Mensal")
lam <- BoxCox.lambda(tsMensalTrain)
lam
## [1] -0.008437
tsMensalBoxCox <- BoxCox(x = tsMensalTrain, lam)
plot(tsMensalBoxCox, col="red")
Calcula-se o fator \(\lambda\) para a transformação e utiliza-se um modelo Exponential Smoothing aditivo (pois os fatores multiplicativos devem ser minimizados pela transformação Box-Cox).
lam <- BoxCox.lambda(tsMensalTrain)
etsMensal_boxcox <- ets(tsMensalTrain, additive=TRUE, lambda=lam)
fcastMensal1yr_boxcox <- forecast(etsMensal_boxcox, h = mesesAPrever)
plot(fcastMensal1yr_boxcox, include=60)
lines(tsMensalTest1yr, col="red")
Acurácia:
accuracy(fcastMensal1yr_boxcox, tsMensalTest1yr)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set 605.1 8996 6818 0.1216 1.456 0.1838 0.008943 NA
## Test set 20121.4 27805 23231 2.6408 3.099 0.6263 0.788149 0.6388
mesesAPrever <- 12
Dados com sazonalidade removida:
tsMensal.stl <- stl(tsMensalTrain[,1], s.window=12)
# Seasonally adjusted data constructed by removing the seasonal component.
plot(seasadj(tsMensal.stl))
stlf_model <- stlf(tsMensalTrain[,1])
stlf_fcast <- forecast(stlf_model, h=mesesAPrever)
plot(stlf_fcast, xlab="Anos", ylab="Fluxo Mensal")
lines(tsMensalTest1yr, col="red")
Acurácia:
accuracy(stlf_fcast, tsMensalTest1yr)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set -120.1 7919 6023 -0.02844 1.297 0.1624 0.02359 NA
## Test set 11044.4 22662 18177 1.36140 2.421 0.4901 0.54398 0.5458
Hyndman, Rob J. 2014. “Forecasting Time Series Using R.” http://robjhyndman.com/talks/MelbourneRUG.pdf.
Hyndman, Rob J, and George Athanasopoulos. 2013. Forecasting: principles and Practice. OTexts.
Hyndman, Rob, Anne B. Koehler, J. Keith Ord, and Ralph D. Snyder. 2008. Forecasting with Exponential Smoothing: The State Space Approach (Springer Series in Statistics). 2008th ed. Springer. http://amazon.com/o/ASIN/3540719164/.
Leek, Jeffrey. 2014. “Practical Machine Learning - Lecture 27 - Forecasting.” Johns Hopkins Bloomberg School of Public Health. https://d396qusza40orc.cloudfront.net/predmachlearn/027forecasting.pdf.